50th AIAA Aerospace Sciences Meeting, January 9-12, 2012, Nashville, Tennessee 


Algebraic Turbulence-Chemistry Interaction Model 


A. T. Norris * 

NASA Langley Research Center, Hampton, Virginia, 23681-2199, U.S.A 


The results of a series of Perfectly Stirred Reactor (PSR) and Partially Stirred Reactor 
(PaSR) simulations are compared to each other over a wide range of operating conditions. 
It is found that the PaSR results can be simulated by a PSR solution with just an adjusted 
chemical reaction rate. A simple expression has been developed that gives the required 
change in reaction rate for a PSR solution to simulate the PaSR results. This expression is 
the basis of a simple turbulence-chemistry interaction model. The interaction model that 
has been developed is intended for use with simple one-step global reaction mechanisms 
and for steady-state flow simulations. Due to the simplicity of the model there is very little 
additional computational cost in adding it to existing CFD codes. 


I. Introduction 


Due to the non-linear nature of the governing equations for turbulent reacting flow, it is well known 
that the expression for the mean species source terms due to chemical reaction are not being calculated 
exactly. Specifically, using the mean species concentrations and mean temperature to calculate the mean 
species source term is not accurate. Indeed it has been shown that the differences in the mean reaction rate 
associated with this approximation can be several orders of magnitude different from the exact solution.^ 
To this end, several models and numerical approaches have been developed in order to resolve this issue. 
Models include the Magnussen model, ^ the Direct Quadrature Moment Method (DQMOM),^ the Assumed 
PDF,^ the Lagrangian PDF^ and the Linear Eddy Model. ^ All these models attempt to model the effect 
that turbulence has on chemical reaction with varying degrees of fidelity. Unfortunately, as the modeling 
becomes more sophisticated, there is a corresponding increase in the computational cost. Of all the models 
listed above, the one with the least computational expense is the Magnussen model, and as a result, it is 
very popular. However the simplicity of the model does result in it lacking several desirable features. 

Let us consider a chemical reaction defined by the mass fraction of Fuel Ij, Oxidizer Yq and Products 
Yp. The fuel reaction rate for the Magnussen model is given by: 


d^l 

dt 


= Min < 


AYfIC 
AYolUrs 
B Yp/tt{l + Tg) 


( 1 ) 


where A and B are constants, Vg is the stoichiometric ratio of fuel to oxidizer and U is the turbulent time 
scale. It is apparent that in this form, the effect of the turbulence on chemistry is modeled by simple making 
the reaction rate inversely proportional to the turbulent time scale and thus as the turbulent time scale 
changes, so does the reaction rate. 
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However real chemical reactions have an upper limit to the reaction rate, which is a not accounted for in 
this model. A simple modification to fix this is to use a global mechanism to limit the reaction rate, such as 
those developed by Westbrook and Dryer. ^ Thus the reaction rate would be give by 


d-Yf ^ I dYf/dt (Magnussen) 

dt \ CYfY^exp(-E/RT) 


( 2 ) 


where C, a, 6, E and R are the global rate constants and T is the temperature. However this simple truncation 
is a somewhat ad-hoc approach and a more rigorous approach would be desirable. 

The goal of this paper is to investigate whether such a model can be developed that has the numerical 
simplicity of the Magnussen model yet have some of the desirable physical effects of the more complex 
methods. The objective is not to develop a model that works for complex finite-rate chemistry problems, 
but rather one that gives a good approximation of reacting flow for a minimum level of computational effort. 
A Magnussen- like reaction rate model, or a simple one-step global mechanism would be examples of the 
target chemical models. 

As an approach to obtaining this improved simple model, consider a single computational cell in a finite- 
volume CFD solution with chemical reaction (figure 1) where rrij is the mass of species j contained in the 






rrij , Sj 








Figure 1. Schematic of single computational cell in a finite volume CFD solution. 


cell and Sj is the rate of creation of species j due to chemical reaction. Also shown is the mass flow rate 
of species j entering the cell, and the mass flow rate leaving the cell, The change in the mass 

contained in the cell as a function of time is given by 

(3) 

For a steady state situation, which would correspond to a steady state CFD solution, the left hand side 
of equation 3 is zero, and it is apparent that there is a strong resemblance between this single cell and a 
Perfectly Stirred Reactor, (PSR). The mass contained in the cell, divided by the inflow mass flux gives a 
residence time, R and the composition of the cell can be plotted simply as a function of residence time. 

The results of a typical PSR calculation are shown in figure 2. As the mass flux into the reactor increases, 
the residence time decreases and the mass fraction of products of combustion start to drop off. When the 
inflow becomes sufficiently large, the reaction is unable to be sustained, and the reactor is said to have blown 
out. For a given volume then, the reactor is parameterized by just the residence time and reaction rate. 
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Figure 2. Results of a typical PSR simulation showing the progress variable plotted as a function of the reactor 
residence time tr • 


Likewise, consider a single cell in an Eulerian Probability Density Function (PDF) simulation shown in 
figure 3. 



jrout 


Figure 3. Schematic of a single computational cell in an Eulerian Probability Density Function simulation. 


The composition in the cell is represented by a probability distribution of different values of scalar mass 
fractions, which evolves due to molecular mixing, M, chemical reaction, S and the inflow and outflow 
of distributions of species j, given by and respectively. The evolution of the distribution of all the 
scalar quantities in the cell, as a function of time is given by 

= E i (4) 

j=i ^ 

where the represents the scalar space of species j and the number of scalars is given by Ng. In a similar 
way as with the PSR, if the change of the distribution of scalars, with time is zero, corresponding to a 
steady state solution, the computational cell can be considered to behave the same as a Partially Stirred 

3 of 15 


American Institute of Aeronautics and Astronautics 


Reactor, (PaSR). Like the PSR, the PaSR is parameterized by a residence time scale and a reaction time 
scale, but also has a mixing time scale. 

While equation 3 and equation 4 look very different they are in fact very similar. If the distribution T 
is set to a dirac delta function, and molecular mixing is not considered, equation 3 is recovered for each 
species j. That is, a PSR can be thought of as a PaSR but without the molecular mixing term. Thus 
by comparing the results of a PSR simulation with the mean PaSR quantities, the effect of the interaction 
between chemistry and turbulence can be evaluated in a consistent manner, and the possibility of a simple 
relationship between the results may be possible. 


II. Discussion 

To investigate the relationship between a PSR and a PaSR, the first thing to do is to develop a suitable 
set of reactor models, and then test them over a range of conditions that might be encountered in a reacting 
flow. In this section the model used for the reactor and the chemistry model are described. The results of a 
comprehensive series of simulations that compare the PSR and PaSR results are also described. 


A. Reactor Model 

In order to compare the results of a PSR and PaSR model, it is important to eliminate any discrepancies 
that may be caused by different numerical techniques. To this end, both the PaSR and the PSR were 
modeled using a particle-based stochastic model. The stochastic method uses an ensemble of particles A/^, 
to represent the scalar composition within a reactor. During a time step, a certain number of particles enter 
the reactor, and an equal number leave. The particles are then mixed via a molecular mixing model 
and the resulting composition is then reacted and the process repeats. 

For this study, three molecular mixing models were used: A Perfect Mixing Model (PMM), the Interaction 
and Exchange with the Mean (lEM) modeff and the Modified Curl Model. ^ 

The PMM takes all the particles in the reactor, and sets their scalar values to the corresponding ensemble 
mean value. Eor example, for a given time step the value of particle (j^ changes as follows: 

= 0(t) (5) 


where (j) is the ensemble mean value. The perfect mixing model, as the name suggests is used to model the 
PSR. 

Eor the PaSR, two models are used. The simpler model is the lEM scheme, which assumes the scalar 
value of the particles will approach the scalar ensemble mean value at an exponential rate governed by the 
turbulent time scale, The change in a particle for a time step of St is given by: 


(j)P(t + 5t) = {(j)P(t) - 0)exp(-ic0^) 

(6) 

where is the mixing constant, with a value commonly set to unity. 

The other PaSR molecular mixing model used is the modified Curl model. This is 
model based on the original Curl model. In this model, a pair of particles, and 

from the ensemble, and their compositions are changed to 

a particle interaction 
are randomly selected 

+ 5t) = 4P{t) + - (f>P{t)) 

(7) 

+ 5t) = + ]^r{4P{t) - 

(8) 


where r is a random number with a uniform distribution in the interval (0,1). In this case however the time 
step St is a function of the turbulent time scale and number of particles in the ensemble, and is given by: 


St = 


U 

3W 


(9) 


If more than one pair of particles are selected at once, the St can be scaled appropriately. 
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In the operation of the PSR and the PaSR model, the particles are initialized with fully mixed and 
reacted particles. The particles that are to be input into the reactor each turn are selected with either a 
mean value or from a specified scalar distribution. 

The overall time step used is based on the number of particles that are selected to be input, however the 
mixing and reaction steps are often run at a smaller time step for several sub iterations to ensure certain 
numerical constraints, such as the time step required for the Curl mixing model. For each run, the code keeps 
iterating until the mean scalar quantities in the reactor have converged to a statistically constant value. The 
mean scalar quantities are evaluated using a time- averaging technique to reduce the amount of scatter in 
the results. 

During the simulations reported in this paper, the results for the modified Curl model are the ones 
predominately reported. However the results for the lEM model have also been performed, and show very 
similar results to those of the modified Curl model. 

B. Reaction Model 

The desired reaction model for this study would be a global one-step reaction, where the species are repre- 
sented by a minimum number of scalar variables. For example the Magnussen model can be represented by 
the mass fraction of fuel, oxidizer and products. This can also be reformulated as simply a mixture fraction 
^ and a progress variable Yp. 

For the purposes of this study, a simple polynomial representation for chemical reaction was developed 
based on ^ and Yp. The motivation for this is to have a mechanism that has the same form as a typical 
one-step global mechanism, yet has an algebraic integral solution that will allow a very efficient solution for 
the PSR and PaSR code. This will afford the opportunity to perform many simulations with a minimum of 
computational effort. The polynomial rate equation chosen is given as 

dY 

-^=AY^{Y,{0-Y,) (10) 

where A is the reaction rate constant, and Yp(^) is the maximum value of the progress variable as a function 
of the mixture fraction, In figure 4 the reaction rate of this chemistry model is compared to scaled versions 
of the Magnussen model and a typical one-step global reaction. 



Figure 4. The progress variable reaction rate dYp/dt plotted as a function of the progress variable Yp for three different 
reaction mechanisms . 
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The rate expression in equation 10 can be integrated to obtain an expression for the progress variable as 
a function of time t, 


1 

AYpHO 



Yp 
- Yp 




( 11 ) 


which is solved via a bisection method. For this study, the value of ^4 = 20 was chosen to give a residence 
time for the reactors of between 0.1 and 10. 

The inverse of the maximum reaction rate for a given value of mixture fraction can be used as a charac- 
teristic time for the reaction. Thus the reaction time scale, tx for equation 10 is given by. 


t 


X 


16 


(12) 


C. Numerical Validity 

To establish the numerical validity of the PSR and PaSR solutions it is important to investigate the numerical 
aspects of the solution process. To this end the numerical parameters controlling the solution were varied to 
obtain a range of operability for the reaction code. 

The first parameter varied was the Particle Replacement Ratio (PRR) number, defined in this effort as 
the number of particles entering the reactor per time step relative to the number of particles contained by 
the reactor. Operating parameters chosen for this were Np = 1000 particles and a turbulence time scale of 
tf =0.2. The results for the modified Curl model simulation are shown in figure 5, where it is apparent that 
the effect of the PRR number on the solution is largest at a residence time of = 1. Similar results were 
also obtained for the PSR and RTM simulations. 



Figure 5. The progress variable plotted as a function of the reactor residence time tr showing the effect of varying 
the particle replacement ratio (PUR). Curl mixing model used. 


Plotting the progress variable against the PRR number for all three mixing models, figure 6, shows that 
the solutions for all models are numerically insensitive for a PRR of less that 0.02, so this value will be used 
in all further simulations. It is interesting to note that the Curl and the RTM model do not converge to 
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Figure 6. The progress variable plotted as a function of the Particle Replacement Ratio (PRR). Results shown for 
the PaSR reactor simulations using the modified Curl model and the lEM mixing model and compared to the PSR 
model. 


the same value despite being expected to provide the same level of mixing. The PSR simulation though was 
expected to converge to a different value as it assumes total mixing. 

The number of particles required for an accurate solution was also investigated. Using the same initial 
conditions as above, except with the PRR number set to 0.02, and varying the total number of particles 
gives a result for the Curl model shown in figure 7. It is apparent from this plot that using N = 100 particles 
gives almost identical results to cases using N = 1000 particles. Thus solving cases with N = 100 particles 
was deemed adequate for all further simulations. Similar results were obtained for the RTM mixing model 
and the PSR. 

For the Curl model the number of particles from the ensemble that are selected for mixing is another 
numerical parameter. Tests were performed and the results were found to be insensitive to the number of 
mixing particles selected. For all further simulations, a maximum of 20% of the ensemble particles were 
selected for mixing each time step. If more particles were required, the mixing and reaction were looped over 
multiple times to ensure the 20% limit was kept to. 

D. Physical Variations 

The first physical parameter varied was the reaction time scale and the PSR results are shown in figure 8. As 
expected, the effect of changing the reaction rate is to simply shift the curve along the residence time axis. 
Doubling the reaction time scale doubles the residence time that the reactor requires to maintain a given 
level of combustion. A similar result was obtained for the PaSR simulation, with the added requirement 
that the ratio of the turbulent and reactive time scales was kept constant. 

The next test was to vary the turbulence time and the PaSR results for the Curl model are shown in 
figure 9. It can be seen in this plot that as the turbulence time scale decreases the PaSR result approaches 
the PSR result. Thus there is an limit to which turbulence will enhance the chemical reaction in a PaSR, 
which is an effect that the basic Magnussen model does not account for. In comparing the curves produced 
by different turbulence levels in figure 9 it can be seen that the shapes are fairly similar to that of the PSR 
result, just moved to the right as the turbulence time scale increases. Coupling this with the result observed 
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tr 


Figure 7. Effect of changing the number of particles Np used to represent the PaSR. 
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Figure 8. Results of PSR simulations with different values of reaction rate. 
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Figure 9. Results of Curl model PaSR simulations for different turbulence time scales. PSR result shown for reference. 


in figure 8 where the PSR curve moves proportional to the reaction rate, a scheme based on changing the 
reaction rate of a PSR to simulate the effects of turbulence in a PaSR should be possible. The expression 
for this model is given by: 

^ = (1 + A it A.) 

where Ao is the original reaction rate constant, A is the new modeled value and Df is a constant which has 
a fitted value of Dt = 2.1 for the kinetic mechanism used in this study. 

The reactor results produced by this expression are shown plotted against the PaSR results for a range 
of turbulence time scales in figure 10. It can be seen that this model does provide a reasonably good 
approximation of the PaSR results. It certainly matches the trend of increased turbulence driving the PaSR 
solution to approach the PSR result. It should be noted that the agreement for lower values of Yp is not so 
good, but as simple chemical mechanisms have considerable uncertainty in this region it is not considered a 
significant flaw. 

Thus for a given set of conditions for a PaSR, the same result can be obtained by a PSR by simply ma- 
nipulating the reaction time scale. This result is the basis of the proposed turbulence-chemistry interaction 
model described in this paper. However in a combustion simulation, the composition entering a computa- 
tional cell will not be a fixed value, and so it is important to investigate the effects of varying the inflow 
distributions of scalar quantities. 

To this end, a simple test was performed that varied the mean progress variable inflow value, Yp^. Sample 
PaSR results of Yp^ =0.2 and Yp^ =0.6 are shown in figure II. What is seen here is that as Yp^ approaches 
the fully burnt value, the effect of turbulent mixing becomes smaller. Thus the model can be modified by 
the mean inflow Yp^ as shown: 


A = 


(1 + A(i - ry/y it A.) 


(14) 


where dt is a constant that is found to fit the data with dt = 3.6. Applying this correction gives results 
shown in figure 12 where it can be seen that a reasonably good agreement has been obtained between the 
PaSR results and the model predictions. 
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Figure 10. Curl model PaSR simulation results for different turbulent time scales with modified PSR results overlayed. 
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Figure 11. Effect of varying mean inflow progress variable value on Curl model PaSR. Lower lines represent a Yp = 0.2 
and upper lines are Yp = 0.6. 
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Figure 12. Results of PSR model with variable inflow composition shown overlayed on Curl PaSR model. Lower lines 
represent a. Yp = 0.2 and upper lines are Yp = 0.6. 


In all the tests so far, the inflow scalar quantities values have just been set to a mean value. However in full 
PDF simulations, the value of the inflow scalar quantities do not just have a single value, but are represented 
by a distribution of values. Because the main goal of this effort was to develop a simple, computational 
efficient model, any effect of a distribution of inflow scalars would ideally be small so they can be ignored. 

To test this, a (3 distribution of was used to generate the inlet quantities. With a range of mean 
progress variable values and standard deviations, ayp = 0.02, 0.1, 0.2 and 0.3 to ensure a wide range of 
possible input scenarios. A sample of the PDFs of the distributions are shown plotted in figure 13 for a mean 
normalized progress variable of 0.2. A sample of the PaSR results are shown in figure 14 where the inflow 
conditions are = 0.2 for four different standard deviations and shown for two different turbulent time 
scales. It can be seen that even for large variations in the distribution, the results are relatively insensitive 
to the precise distribution and so this effect can be ignored. 

So far in this study, the effects of varying the mixture fraction have not been addressed, as this is a 
conserved scalar. However there can be an effect on the mean progress variable when there is a distribution 
of To test this, a relationship between the maximum progress variable Ypm and the mixture fraction ^ 
was chosen and is shown plotted in figure 15. By testing at the stoichiometric value = 0.2 it is expected 
that the PaSR results will be the most effected, as all non-stoichiometric values of mixture fraction will 
exhibit lower reaction rates than the mean value of ^ and so drag the results lower. A beta distribution 
of ^ values was chosen to represent possible distributions, with a mean ^ = 0.2 and a standard deviation 
of (7^ = 0.02, 0.05, 0.10 and 0.20 are shown in figure 16. These distributions were used with a variety 
of different combinations of progress variable mean values and distributions and a sample result for two 
different turbulence time scales is shown in figure 17. While there is some variation due to the change in ^ 
distribution, the differences are sufficiently small that they can be ignored for this model. 

What has been shown is that a simple modification (shown in equation 14) to a global reaction rate can 
provide some basic features of a turbulence-chemistry interaction model. For a given global reaction, the 
two constants Bf and bt may have to be calibrated, but the values given here should be a reasonable start 
for most forms of reactions. In addition, due to the simple form of equation 14 the additional computational 
cost caused by this model should be negligable. 
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Figure 13. Selection of the distributions of progress variables used in PaSR simulations with mean of 0.2. B_sig refers 
to the standard deviation of the distribution. 



tr 


Figure 14. Variation of PaSR results due to change of inflow progress variable distributions. Square symbols are for 
tt = 1 and circles are for tt = 10. B_sig refers to the standard deviation of the distribution. 
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Figure 15. Maximum progress variable plotted against mixture fraction. 





Figure 16. Distributions of mixture fraction use in PaSR distribution. B_sig refers to the standard deviation of the 
distribution. 
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Figure 17. Variation of PaSR results due to change of inflow mixture fraction distributions. Square symbols are for 
tf = 1 and circles are for tt = 10. B_sig refers to the standard deviation of the distribution. 


III. Conclusions 

By comparing the mean scalar values obtained from a Perfectly Stirred Reactor (PSR) with those of a 
Partially Stirred Reactor (PaSR) a simple algebraic expression has been obtained that allows some of the 
features of the interaction between chemistry and turbulence to be modeled. 

First a particle-based stochastic method was used to simulate both the PSR and PaSR. Tests were then 
performed to ensure the numerical accuracy of the scheme. While the results were relatively insensitive to 
the total number of particles used to represent the reactor, the number of particles flowing into the reactor 
per time step had to be limited to a small percentage of the total number to ensure accuracy. A parametric 
study of the physical parameters governing the performance of the reactors was then performed. It was 
found that the PaSR results were most sensitive to the variation of the mean inflow composition and the 
turbulent time scale, while the effect of varying the distribution of inflow quantities was small. As a result, 
a simple model has been proposed that models the effect of turbulence and mean inflow composition for a 
reacting flow. 

Implementation of the model requires only the modification of the reaction rate term and the evaluation 
of a reaction rate time scale, defined as the inverse of the peak reaction rate. The model requires two 
constants, that may need to be tuned for different reaction rate expressions, however these values can be set 
via a PaSR simulation and should be constant for all uses of that kinetic mechanism. It should be noted 
that this expression is intended for use only with simple one-step reaction schemes and not with complex 
multi-step mechanisms. 
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